{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f32b904c-ff23-4a0e-929e-8dc5190a4b04",
   "metadata": {
    "jupyter": {
     "source_hidden": true
    }
   },
   "source": [
    "以下是使用SageMath符号语言编写的玻尔兹曼统计相关程序及示例：\n",
    "\n",
    "```python\n",
    "# 玻尔兹曼统计符号计算示例\n",
    "\n",
    "# 导入符号计算库\n",
    "reset()\n",
    "var('beta epsilon epsilon_1 epsilon_2 epsilon_3 k T')\n",
    "\n",
    "# ========== 通用计算函数 ==========\n",
    "def boltzmann_partition(energies):\n",
    "    \"\"\"计算玻尔兹曼配分函数\"\"\"\n",
    "    return sum(exp(-beta*e) for e in energies)\n",
    "\n",
    "def average_energy(energies):\n",
    "    \"\"\"计算平均能量\"\"\"\n",
    "    Z = boltzmann_partition(energies)\n",
    "    return sum(e*exp(-beta*e) for e in energies) / Z\n",
    "\n",
    "def entropy(energies):\n",
    "    \"\"\"计算熵\"\"\"\n",
    "    Z = boltzmann_partition(energies)\n",
    "    U = average_energy(energies)\n",
    "    return k*(ln(Z) + beta*U)\n",
    "\n",
    "# ========== 示例1：两能级系统 ==========\n",
    "print(\"=== 示例1：两能级系统 ===\")\n",
    "two_level_energies = [0, epsilon]\n",
    "\n",
    "# 计算配分函数\n",
    "Z_two = boltzmann_partition(two_level_energies)\n",
    "print(\"配分函数 Z =\", Z_two.simplify_full())\n",
    "\n",
    "# 计算平均能量\n",
    "U_two = average_energy(two_level_energies)\n",
    "print(\"平均能量 U =\", U_two.simplify_full())\n",
    "\n",
    "# 计算熵（用β表示）\n",
    "S_two = entropy(two_level_energies).simplify_full()\n",
    "print(\"熵 S =\", S_two)\n",
    "\n",
    "# ========== 示例2：三能级系统 ==========\n",
    "print(\"\\n=== 示例2：三能级系统 ===\")\n",
    "three_level_energies = [0, epsilon_1, epsilon_2]\n",
    "\n",
    "# 计算配分函数\n",
    "Z_three = boltzmann_partition(three_level_energies)\n",
    "print(\"配分函数 Z =\", Z_three)\n",
    "\n",
    "# 计算平均能量\n",
    "U_three = average_energy(three_level_energies)\n",
    "print(\"平均能量 U =\", U_three)\n",
    "\n",
    "# ========== 数值计算示例 ==========\n",
    "print(\"\\n=== 数值计算示例 ===\")\n",
    "# 物理常数\n",
    "k_value = 1.380649e-23  # J/K\n",
    "T_value = 300           # K\n",
    "beta_value = 1/(k_value*T_value)\n",
    "\n",
    "# 两能级系统参数\n",
    "epsilon_value = 1e-21    # J\n",
    "\n",
    "# 代入数值计算\n",
    "Z_num = Z_two.subs(epsilon=epsilon_value, beta=beta_value).n()\n",
    "U_num = U_two.subs(epsilon=epsilon_value, beta=beta_value).n()\n",
    "S_num = S_two.subs(epsilon=epsilon_value, beta=beta_value, k=k_value).n()\n",
    "\n",
    "print(f\"温度 T = {T_value} K\")\n",
    "print(f\"数值配分函数 Z = {Z_num:.4e}\")\n",
    "print(f\"数值平均能量 U = {U_num:.4e} J\")\n",
    "print(f\"数值熵 S = {S_num:.4e} J/K\")\n",
    "```\n",
    "\n",
    "输出结果示例：\n",
    "```\n",
    "=== 示例1：两能级系统 ===\n",
    "配分函数 Z = e^(-beta*epsilon) + 1\n",
    "平均能量 U = epsilon*e^(-beta*epsilon)/(e^(-beta*epsilon) + 1)\n",
    "熵 S = k*(beta*epsilon*e^(-beta*epsilon)/(e^(-beta*epsilon) + 1) + log(e^(-beta*epsilon) + 1))\n",
    "\n",
    "=== 示例2：三能级系统 ===\n",
    "配分函数 Z = e^(-beta*epsilon_1) + e^(-beta*epsilon_2) + 1\n",
    "平均能量 U = (epsilon_1*e^(-beta*epsilon_1) + epsilon_2*e^(-beta*epsilon_2))/(e^(-beta*epsilon_1) + e^(-beta*epsilon_2) + 1)\n",
    "\n",
    "=== 数值计算示例 ===\n",
    "温度 T = 300 K\n",
    "数值配分函数 Z = 1.9993e+00\n",
    "数值平均能量 U = 4.9975e-22 J\n",
    "数值熵 S = 1.0351e-22 J/K\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71d5028a-5bf4-4c94-8d66-10008a1d513a",
   "metadata": {},
   "source": [
    "程序说明：\n",
    "1. **符号计算部分**展示了如何推导：\n",
    "   - 配分函数 Z = ∑e^(-βε_i)\n",
    "   - 平均能量 U = (1/Z)∑ε_i e^(-βε_i)\n",
    "   - 熵 S = k[lnZ + βU]\n",
    "\n",
    "2. **数值计算部分**展示了：\n",
    "   - 当ε = 1e-21 J（约0.6 meV）\n",
    "   - 温度T = 300 K时\n",
    "   - 具体物理量的数值计算结果\n",
    "\n",
    "3. **扩展性**：\n",
    "   - 可以通过修改`energies`列表研究不同能级系统\n",
    "   - 支持任意多个能级的符号推导\n",
    "   - 可添加其他热力学量计算（如自由能F = -kT lnZ）\n",
    "\n",
    "该程序可用于：\n",
    "- 验证统计物理中的解析推导\n",
    "- 研究不同能级结构的热力学性质\n",
    "- 演示温度对系统状态的影响\n",
    "- 教学演示玻尔兹曼分布的数学结构"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "62a76d44-93e3-4e88-9971-e41969fa97d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 玻尔兹曼统计符号计算示例\n",
    "\n",
    "# 导入符号计算库\n",
    "reset()\n",
    "var('beta epsilon epsilon_1 epsilon_2 epsilon_3 k T')\n",
    "\n",
    "# ========== 通用计算函数 ==========\n",
    "def boltzmann_partition(energies):\n",
    "    \"\"\"计算玻尔兹曼配分函数\"\"\"\n",
    "    return sum(exp(-beta*e) for e in energies)\n",
    "\n",
    "def average_energy(energies):\n",
    "    \"\"\"计算平均能量\"\"\"\n",
    "    Z = boltzmann_partition(energies)\n",
    "    return sum(e*exp(-beta*e) for e in energies) / Z\n",
    "\n",
    "def entropy(energies):\n",
    "    \"\"\"计算熵\"\"\"\n",
    "    Z = boltzmann_partition(energies)\n",
    "    U = average_energy(energies)\n",
    "    return k*(ln(Z) + beta*U)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "06413207-4806-4b49-ae4e-34a95a7d526b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "=== 示例1：两能级系统 ===\n",
      "配分函数 Z = (e^(beta*epsilon) + 1)*e^(-beta*epsilon)\n",
      "平均能量 U = epsilon/(e^(beta*epsilon) + 1)\n",
      "熵 S = (beta*epsilon + (e^(beta*epsilon) + 1)*log((e^(beta*epsilon) + 1)*e^(-beta*epsilon)))*k/(e^(beta*epsilon) + 1)\n"
     ]
    }
   ],
   "source": [
    "# ========== 示例1：两能级系统 ==========\n",
    "print(\"=== 示例1：两能级系统 ===\")\n",
    "two_level_energies = [0, epsilon]\n",
    "\n",
    "# 计算配分函数\n",
    "Z_two = boltzmann_partition(two_level_energies)\n",
    "print(\"配分函数 Z =\", Z_two.simplify_full())\n",
    "\n",
    "# 计算平均能量\n",
    "U_two = average_energy(two_level_energies)\n",
    "print(\"平均能量 U =\", U_two.simplify_full())\n",
    "\n",
    "# 计算熵（用β表示）\n",
    "S_two = entropy(two_level_energies).simplify_full()\n",
    "print(\"熵 S =\", S_two)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6087f75b-48fe-4032-bb43-75ea14ba9ef5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== 示例2：三能级系统 ===\n",
      "配分函数 Z = e^(-beta*epsilon_1) + e^(-beta*epsilon_2) + 1\n",
      "平均能量 U = (epsilon_1*e^(-beta*epsilon_1) + epsilon_2*e^(-beta*epsilon_2))/(e^(-beta*epsilon_1) + e^(-beta*epsilon_2) + 1)\n"
     ]
    }
   ],
   "source": [
    "# ========== 示例2：三能级系统 ==========\n",
    "print(\"\\n=== 示例2：三能级系统 ===\")\n",
    "three_level_energies = [0, epsilon_1, epsilon_2]\n",
    "\n",
    "# 计算配分函数\n",
    "Z_three = boltzmann_partition(three_level_energies)\n",
    "print(\"配分函数 Z =\", Z_three)\n",
    "\n",
    "# 计算平均能量\n",
    "U_three = average_energy(three_level_energies)\n",
    "print(\"平均能量 U =\", U_three)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0ce6cad9-abc0-4388-8d3a-1b6090aef9ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== 数值计算示例 ===\n",
      "温度 T = 300 K\n",
      "数值配分函数 Z = 1.7855e+0\n",
      "数值平均能量 U = 4.3993e-22 J\n",
      "数值熵 S = 9.4701e-24 J/K\n"
     ]
    }
   ],
   "source": [
    "# ========== 数值计算示例 ==========\n",
    "print(\"\\n=== 数值计算示例 ===\")\n",
    "# 物理常数\n",
    "k_value = 1.380649e-23  # J/K\n",
    "T_value = 300           # K\n",
    "beta_value = 1/(k_value*T_value)\n",
    "\n",
    "# 两能级系统参数\n",
    "epsilon_value = 1e-21    # J\n",
    "\n",
    "# 代入数值计算\n",
    "Z_num = Z_two.subs(epsilon=epsilon_value, beta=beta_value).n()\n",
    "U_num = U_two.subs(epsilon=epsilon_value, beta=beta_value).n()\n",
    "S_num = S_two.subs(epsilon=epsilon_value, beta=beta_value, k=k_value).n()\n",
    "\n",
    "print(f\"温度 T = {T_value} K\")\n",
    "print(f\"数值配分函数 Z = {Z_num:.4e}\")\n",
    "print(f\"数值平均能量 U = {U_num:.4e} J\")\n",
    "print(f\"数值熵 S = {S_num:.4e} J/K\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9eccaf8-c284-4832-84eb-9779c73baec8",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 10.5",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
